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ABSTRACT 

We propose a new method to conserve the total energy to round-off error in grid- 
based codes for hydrodynamic simulations with self-gravity. A formula for the energy 
flux due to the work done by the the self-gravitational force is given, so the change 
in total energy can be written in conservative form. Numerical experiments with the 
code Athena show that the total energy is indeed conserved with our new algorithm 
and the new algorithm is second order accurate. We have performed a set of tests that 
show the numerical errors in the traditional, non-conservative algorithm can affect the 
dynamics of the system. The new algorithm only requires one extra solution of the 
Poisson equation, as compared to the traditional algorithm which includes self-gravity 
as a source term. If the Poisson solver takes a negligible fraction of the total simulation 
time, such as when FFTs are used, the new algorithm is almost as efficient as the 
original method. This new algorithm is useful in Eulerian hydrodynamic simulations 
with self-gravity, especially when results are sensitive to small energy errors, as for 
radiation pressure dominated flow. 

Subject headings: methods: numerical — gravitation — hydrodynamics 



1. Introduction 



Self-gravity plays an important role in many astrophysical flows. For instance, it is central to 
massive star formation (e.g., McKee & Ostriker 2007; Krumholz et al. 2009), supernova explosions 
(e.g., Nordhaus et al. 2010), planet formation in protoplanetary disks (e.g., Armitage 2011), and 
structure formation in the early universe (e.g., Bertschinger 1998). Numerical modeling of these 
systems requires that self-gravity be implemented correctly in hydrodynamical simulations, espe- 
cially in problems for which energy balance is important. For example, the fate of a collapsing 
cloud depends on the total energy which is the sum of the potential, internal, and kinetic energies. 
For some systems, e.g. those having a relativistic equation of state with adiabatic gamma, 7 ~ 4/3, 
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the total energy can be much smaller than either the potential or internal energies in hydrostatic 
equilibrium. Thus, a small numerical error in the computation of the potential energy can lead 
to a large fractional error in the total energy, potentially causing a bound system to become un- 
bound (e.g., Jiang & Goodman 2011) or vice versa. Energy conservation in a grid-based code 
with self-gravity is usually not guaranteed. The goal of this paper is to propose a new algorithm 
for grid-based codes, which conserves the sum of the internal, kinetic, and potential energies to 
round-off error, allowing accurate simulations of self-gravitating systems to be performed. 

The change in the momentum and the energy due to self-gravity are typically added to the 
equations of hydrodynamics as source terms. In this case, the Euler equations read 

^ + V-(p») = 0, (1) 

^ + V-(pvv + P) = -pVcp, (2) 
r) W 

— + V-[(E + P)v)] = -pv-V<f>, (3) 

V 2 4> = AttG P . (4) 

Here p, v,P are density, velocity, and pressure respectively, and E is the sum of the internal and 
kinetic energies. G is the gravitational constant, and (j) is the gravitational potential, which is 
related to the density via the Poisson equation (4) . We can also define the total energy as 

Etot = E + \pcj). (5) 



E to t is a globally conserved quantity so that 

d_ 

di 



J d 3 X E tot = (6) 



for any isolated self-gravitating system. Expression (6) can be derived by integrating equation (3) 
over all of space and applying the divergence theorem (see e.g., Binney Sz Tremaine 2008, section 
2.1). 

Equations (2) and (3) are written in non-conservative form. Therefore, momentum and energy 
are usually not conserved to roundoff error when they are solved numerically. It is well known that 
if momentum is not conserved for a self-gravitating system, then the shape of the object cannot be 
kept very well when it moves across the simulation box (e.g., Edgar et al. 2005; Tasker et al. 2008). 
In order to conserve momentum, it is now common practice to write the momentum equation in 
conservative form as 

^ + V-(pvv + P + T g ) = 0, (7) 
where the gravitational tensor T g is 



T - J- 



V^V0-i(V<A)-(Vc 



(8) 
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and I is the identity tensor. This equation is mathematically equivalent to the original momentum 
equation (2) with use of equation (4). However equation (8) can be solved in a completely different 
way numerically. Momentum change due to gravitational acceleration is no longer included as a 
source term. Instead, it is included as a flux term via T g , which means the total momentum can be 
conserved to round-off error numerically. This has already been implemented in many grid-based 
codes, such as Athena. 

The improvement we propose in this paper is to also write the energy equation with self-gravity 
(3) in conservative form. Using this form of the equation, we devise a numerical algorithm that 
conserves the total energy with self-gravity to round-off error. 



2. Formula for the energy flux 

We begin by rewriting the energy equation (3) in the following form 

Wti E + \ p<t ) + v ' [{E + p)v + Fg)] = °- (9) 

As long as the flux due to self-gravity F g falls off faster than r -2 at large distances and the system 
is spatially bounded, equation (9) automatically implies that equation (6) is satisfied. We now 
solve for the form of F g . 

Using equation (3), equation (9) can be written as 

Ft{l P<P ) +V-F a = pvVct>. (10) 

Using the continuity equation, this becomes 

V • F g = - V • (pv) cp + pv-Vcfy- -p$. (11) 

Here <fi = d<p/dt is the time rate of change of the potential. Differentiating the Poisson equation 
with respect to time yields 

V 2 = 4vrGp = -4vrGV • (pv) . (12) 
Substituting this in equation (11) gives 

V-F g = V • (pv) cp + pv ■ + ^ (V>V 2 - </)V 2 0) (13) 

= v- pw> + ^ (<W - <w>) 

so one form of the energy flux due to self-gravity is 

F 9 = ^Q ( 0V< ^ " ^ V0 ) + ^ (14) 
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Note that this form for the energy flux is not unique, and as discussed in Appendix A, alternative 
forms are available, which also satisfy equation (13). However, only the divergence of the gravita- 
tional tensor and energy flux are used to evolve the system. The gravitational tensor and energy 
flux themselves are not used directly. Thus, as long as equation (13) is satisfied, different forms of 
the energy flux are mathematically equivalent. We note that a similar formula for energy flux due 
to self-gravity has been proposed by Pen (1998). 

In order to calculate the energy flux, we need to solve equation (12) for cf>. Note that this 
equation has the same form as the Poisson equation, which means it can be solved with the same 
numerical technique (such as FFT for periodic boundary conditions). Details on how to calculate 
the energy flux F g are given in the next section. 



3. Numerical Implementation 

We implement the equations of hydrodynamics with self-gravity in conservative form using 
Athena (e.g., Stone et al. 2008), which is a MHD code that uses an unsplit, higher order Godunov 
scheme. Since the implementation of self-gravity is unchanged by the addition of a magnetic field, 
we focus on the hydrodynamic case here. With self-gravity, the following steps are needed in 
addition to the original algorithm described in Stone et al. (2008). 

First, at time step n we need to compute the gravitational potential <j) n by solving the Pois- 
son equation using the density distribution. For periodic boundary conditions, this can be done 
efficiently using the Fast Fourier Transform (FFT). 

Second, after we get the left and right states at the cell interface for each direction (step 1 
of 6.1 in Stone et al. 2008), we need to add the change due to self-gravity to the left and right 
states at the half time step 5t/2. If primitive variables (density, velocity and pressure) are used for 
the reconstruction, we only need to add the gravitational acceleration — to the left and right 
velocity. No gravitational source terms need to be added to the pressure. 

Third, after we get the fluxes at the interfaces from the Riemann solvers (step 7 in Stone 
et al. 2008), we solve equation (12) for <fi n+1 / 2 , and the same numerical technique as for solving 
the Poisson equation can be used. We just need to replace p on the right hand side of the Poisson 
equation with —V • (pv), where pv is the density flux from the Riemann solvers at each interface. 
The 4> n+l / 2 we get is at the cell centers and is the time-averaged value for time step n. 

Fourth, we update the density from time step nton + 1 and calculate the new potential (f) n+l 
based on the updated density distribution. Then, we use the averaged potential (4> n + (f) n+l )/2 and 
0n+i/2 to ca i cu i a te the gravitational tensor T g and the energy flux F g . The cell-centered potential 
4> is spatially-averaged to get the potential at cell interfaces. 

Finally, we update the momentum and the total energy with T g and F g . At each cell (i, j, k), 
if the cell size is 5x x 5y x 5z, then the momentum change (take pv x for example) due to self-gravity 
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is 



/ T n+1/2 _ T n+1/2 
£ x ^ 9xxi+l/2,j,A; J flxxi-l/2, ; 

^ / T n+1/2 

<Sy V B v*ij+i/'. 



and the change in the total energy 



— (t 71+1,2 

fi z \ 9zxi,j,k+l 



t+1/2 _ j, n+1/2 

,j+l/2,k 9 y*i,j-i/: 



n+1/2 _ j, n+1/2 
zxi,j,k+l/2 9zxi,j,k—l 



l/2,j,k) 

,j-l/2,k) 

+1/2 \ 
zxi,j,k-l/2J - 



is 



rpn+1 T?n , 

^i,j,k ~ Hj,k + 



2 " 2 (p ^K 



^ / F n+1/2 

5x V 9 ^+l/2,; 



re+1/2 p n+1/2 \ 

i+l/2,j,k r 9xi-l/2,j,k) 



_ (p n+l / 2 - F n+1/2 
5y V 9 ^,j+i/2,fc 9 Vi,j-i/2,k) 

n+1/2 w n+1/2 \ 

§z \ r 9zi,j,k+l/2 r 9zi,j,k-l/2) ■ 

Note that the definition of E does not include potential energy. 



(15) 
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4. Numerical Tests 



To show that our algorithm is stable, second order accurate, and conserves the total energy to 
machine precision, we perform several tests in one, two, and three dimensions. We also compare 
the results of the new algorithm with a non-conservative one that has self-gravity added as a source 
term. In all cases, the non-conservative algorithm use the conservative form of the momentum - 
but not energy — equation. In the non-conservative algorithm, the gravitational potential 4> n is first 
calculated based on density at time step n. After the density is updated, the new potential cf) n+1 at 
time step n+l is calculated. Then the energy source term is added as — (pv) n+1 ^ 2 - V (<f> n + <p n+1 ) /2, 
where (pv) n+1 ^ 2 is the average of the left and right sides of the density flux. In this way, the non- 
conservative algorithm can achieve second order accurate as long as energy error due to self-gravity 
is negligible. 



4.1. Jeans Collapse in ID and 2D 

All of our tests in ID and 2D are based on the Jeans instability problem, where self-gravity is 
the dominant force for the dynamics. We initialize the background state to be a self-gravitating, 
uniform density medium and assume the fluid is adiabatic with sound speed c 2 = jP/p. On top 
of this background state, we initialize a small-amplitude (Sp/p = 1CT 6 ) normal mode perturbation 
with wavevector k. The amplitude of the perturbation changes with time as e tu)t ', where u is given 
by the Jeans dispersion relation 



oj 2 = k 2 c 2 s - AnGp. 



(17) 
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The Jeans length is given by Aj = y/irc^/Gp, so that a normal mode perturbation with wavelength 
A < Aj results in a propagating wave, whereas one with A > Aj yields an instability that saturates 
by forming dense clumps. We focus on the unstable case for comparing the conservative and non- 
conservative algorithms. Table 1 lists the details of each of our simulations and shows the relevant 
numerical parameters used. 



The first tests we perform are an accuracy and a convergence test in ID. In Figure 1, the left 
panel shows the time-evolution of the square of the perturbed gravitational potential integrated over 
the simulation domain for simulation A, which is for an unstable mode having A = 2Aj. The points 
show values from the simulations with time measured in units of {kjc s )~ l and distance measured 
in units of Aj. The gravitational potential is thus measured in units of (Xjkjc s ) 2 = (2irc s ) 2 . The 
solid line shows the analytical solution obtained from the dispersion relation (17). It is clear that 
the code accurately captures the exponential growth in the linear phase of the instability through 
to saturation. The right panel of Figure 1 shows convergence of the new algorithm with resolution. 
The initial perturbation has A = Aj/2 and corresponds to a stable perturbation. The points are for 
simulations B-F (the only thing that varies between the simulations is resolution) and the x-axis 
shows the number of cells per Jeans length. The y-axis shows the error in the Li norm between 
the numerical and analytical solutions for $ after one period (t = 2ir/u). The error in the £2 norm 
can be represented as 



where 3>j is the value of the exact solution for the potential at point i, <I>j is the value of the 
numerical solution, and the summation runs over all points in the simulation domain. For reference, 
we also plot a line whose slope indicates second order convergence, and which confirms that the 
new algorithm does indeed converge at second order. 

To test the conservation properties of the new algorithm in comparison with the old one, we 
compute the fractional energy error at timestep n to be \E n — E°\/E°, where 



Indeed, our new algorithm conserves the total energy to round-off error. However, in ID, the 
original algorithm also conserves the total energy to very high precision. Thus, we need to go to 
a higher number of dimensions to demonstrate the superior energy conservation properties of the 
new algorithm. 



4.1.1. ID tests 




(18) 




(19) 
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sim 


dim 


alg 


N 


A 


Aj 


M 


7 


A 


1 


c 


256 


256 


128 





5/3 


B 


1 


c 


16 


16 


32 





5/3 


B2 


2 


c 


23x23 


16 


32 





5/3 


C 


1 


c 


32 


32 


64 





5/3 


C2 


2 


c 


45x45 


32 


64 





5/3 


D 


1 


c 


64 


64 


128 





5/3 


D2 


2 


c 


91x91 


64 


128 





5/3 


E 


1 


c 


128 


128 


256 





5/3 


E2 


2 


c 


181x181 


128 


256 





5/3 


F 


1 


c 


256 


256 


512 





5/3 


F2 


2 


c 


362x362 


256 


512 





5/3 


G 


1 


c 


512 


512 


1024 





5/3 


H 


2 


nc 


22x22 


16 


8 





5/3 


I 


2 


nc 


45x45 


32 


16 





5/3 


J 


2 


nc 


91x91 


64 


32 





5/3 


Jc 


2 


c 


91x91 


64 


32 





5/3 


K 


2 


nc 


181x181 


128 


64 





5/3 


L 


2 


nc 


362x362 


256 


128 





5/3 


M 


2 


nc 


91x91 


64 


32 





4/3 


N 


2 


nc 


181x181 


128 


64 





4/3 





2 


nc 


362x362 


256 


128 





4/3 


P 


2 


c 


91x91 


64 


32 


10 


5/3 



Table 1: From left to right, the columns designate the simulation label, the number of dimensions, 
the algorithm used (c for conservative, nc for nonconservative) , the number of cells (N x x N y x N z ), 
the wavelength of the initial perturbation measured in cells (all cells have the same size), the Jeans 
wavelength measured in cells, the Mach number of the background fluid, and the adiabatic index 
of the gas. The Courant number for all simulations is 0.8, and the HLLC Riemann solver is used 
for all simulations. 
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Fig. 1. — Left: Comparison of the analytical and numerical solutions for simulation A. The an- 
alytical solution for the linear stage of growth is shown by the line, and the numerical solution 
is shown with points. At late times, the numerical solution deviates from the linear analytical 
solution because the instability becomes non-linear and saturates. Right: Convergence study for 
the conservative algorithm. The points correspond to simulations B — F, and a line with slope 
indicating second order convergence is plotted for reference. 

4-1-2. 2D linear wave tests 

To see how our new algorithm and the non-conservative algorithm behave in 2D, we again 
initialize an unstable eigenmode with A = 2Aj, just as in the ID case. Now, however, we choose 
k x /ky = 1 so that the wavefronts are not aligned with the grid. 

We perform two kinds of tests for the 2D case, which are analogous to the tests in the ID case. 
The first is a convergence test of the conservative and non-conservative algorithms, for a stable 
mode (A < Aj). The second, is a test of the energy conservation for the two algorithms for an 
unstable mode (A > Aj). 

Figure 2 shows the results of the convergence test for the conservative algorithm (simulations 
B2-F2). On a log-log graph, we plot the £2 norm error in p as a function of resolution after one 
oscillation period for a mode having A = Aj/2, and an initial amplitude of 5p/p = 10 -6 . The points 
denote the simulation results, and the slope of the solid line shows second order convergence. It is 
clear that the conservative algorithm does indeed converge at second order. We have tested that 
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the non-conservative algorithm also converges at second order for this the stable mode. Because 
the results are essentially identical to the conservative algorithm, we do not plot them. 




number of gridpoints per Jeans length 



Fig. 2. — Convergence study for the conservative algorithm in 2D. The points correspond to simu- 
lations B2 — F2, and a line with slope indicating second order convergence is plotted for reference. 
The non-conservative algorithm gives almost the same result for this stable propagation linear 
mode. 

We now describe tests of energy conservation, using an unstable normal mode. In 2D, in order 
to make the normal mode collapse to a filament as opposed to a sheet, we increase the density 
perturbation in the center of the simulation box by 2% relative to the edges. We find that the 
energy conservation properties of the non-conservative algorithm are drastically degraded relative 
to the ID case. Table 2 shows the fractional energy error for the nonconservative 2D, 7 = 5/3 
simulations (H — L), and the nonconservative 2D, 7 = 4/3 simulations (M — O). It is clear 
from the table, that in 2D the energy error for the non-conservative algorithm can be of order the 
total energy, although the error does decrease with increasing resolution. On the other hand, the 
conservative algorithm still conserves the total energy to machine precision for the same set of tests. 

Note that the initial total energy (ignoring the energy of the perturbation) is roughly E = 
P/(7 — 1). Thus, the initial total energy for the 7 = 4/3 case is twice that for the 7 = 5/3 case. 
This means that although the fractional energy error for the 7 = 5/3 case is roughly double that 
of the 7 = 4/3 case, the absolute value of the accumulated energy error is roughly the same. We 
have also made sure that the total energy is conserved for the conservative algorithm if we advect 
the background medium with respect to the grid at Mach 10 (simulation P). 
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Label Error 



Label Error 



H .65 

I .32 

J .18 

K .09 

L .05 



M 



N 



O 



.09 
.05 



.03 



Table 2: The energy error for the non-conservative simulations in 2D as a function of resolution 
(increasing resolution going down). The error shown is the fractional error in the total energy 
between the beginning of the simulation and when the perturbation has virialized. The simulations 
H, I, J, K, L are done with adiabatic index 7 = 5/3 while the simulations M, N, O are done with 
adiabatic index 7 = 4/3. Note that in the case 7 = 4/3, the initial total energy is twice that of the 
case 7 = 5/3 (neglecting the energy of the perturbation). 

The main cause of the discrepancy between the ID and the 2D case is likely to be twofold. 
First, the wavevector is no longer aligned with the grid, so the 2D case is less symmetric than the 
ID case. Second, the collapse proceeds further in the 2D case: to a filament in 2D as opposed 
to a sheet in ID. As self-gravitational energy is a global quantity, the energy error caused by the 
source terms of self-gravity has different dependence on time step and grid size compared with the 
normal local truncation error. This energy error will be at a maximum when the gravitational 
potential changes most rapidly When the total numerical error is dominated by the gravitational 
energy error for the non-conservative algorithm, Table 2 shows that the error scales roughly as 
0(1/N). However, when the energy error from self-gravitational source terms is much smaller than 
the usual truncation error as in the case of Figure 2, the non-conservative algorithm shows second 
order convergence with resolution. 

Fig. 3 shows the time-evolution of the energy error for simulation H. Overplotted is the kinetic 
energy in the simulation, normalized to the maximum kinetic energy. It is evident that most of the 
the error is accumulated in a short interval around the time when the kinetic energy, and hence the 
infall velocity are large. This is also when the gravitational potential changes most rapidly, making 
the energy error large for the non-conservative algorithm. 

For reference, Fig. 4 shows three snapshots of the density from the 2D simulation which has 
a resolution of 128 cells per Jeans length. The first snapshot is taken during the linear phase of 
exponential growth, the second during the non- linear phase when matter is collapsing to a sheet, 
the third during the phase when matter is collapsing to a filament, and the fourth during the final, 
virialized state. 

In order to see the possible effect on the dynamics due to the energy error, we compare 
the density distribution from simulations J and Jc. The two simulations have exactly the same 
parameters but simulation J uses the non-conservative algorithm, whereas simulation Jc uses the 
conservative one. Figure 5 shows that the density distributions are different depending on whether 
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Fig. 3. — The solid line shows the fractional error in the total energy as a function of time for 
simulation H. The dashed line shows the kinetic energy as a fraction of the maximum kinetic 
energy in the simulation. It is clear that the error accumulates when the kinetic energy is highest. 
The plateau between t = 10 and t = 20 corresponds to the phase in which the collapsing object 
has already finished collapsing to a sheet, but is just starting to collapse to a filament. 

the total energy is conserved or not. When the total energy is conserved (right panel of Figure 5), 
the formed clump is symmetric, as it should be, because it is formed from a symmetric filament. 
However, for the non-conservative algorithm, the symmetry is lost (left panel of Figure 5). Moreover, 
rotation is generated with the non-conservative algorithm, which shouldn't be the case, since the 
initial angular momentum is zero. This is likely because the energy error of the non-conservative 
algorithm changes the distribution of gas pressure, which causes the rotation of the fluid. 



The test that shows the difference between the conservative and non-conservative algorithms 
most clearly in 3D is the collapse of a polytropic sphere. This test is also relevant to many 
astrophysical systems where self-gravity is important, such as star formation and planet formation 
in protoplanetary disks. In 3D, the adiabatic index 7 = 4/3 is a critical case when the total energy 
of a hydrostatic self-gravitating sphere becomes zero. The dynamics of the sphere will be sensitive 
to the energy error made in the numerical solution, which may change the sign of the total energy in 
the solution. The 7 = 4/3 case is relevant for massive stars when gravity is balanced by a radiation 
pressure gradient (e.g., Jiang & Goodman 2011). 



4.2. Collapse of a Polytropic Sphere in 3D 



We initialize a polytropic sphere with polytropic index n = 3 in a periodic domain, so the 




Fig. 4. — Snapshots of the density and the velocity field from simulation L. The scale of the density 
for each figure is given by the corresponding color bar, and the initial background density is p = 1. 
The magnitude scale of the velocity vectors is arbitrary, and varies from figure to figure. The upper 
left panel is taken at t = 9.2 (A;jc s ) _1 , during the linear stage of the collapse. The upper right panel 
shows the simulation at t = 18 (/cjc s ) _1 , in the non-linear regime when the density distribution 
has collapsed to a sheet. The lower left panel shows the collapse of the sheet to a filament at 
t = 21 (kjCg)" 1 . The lower right panel shows the final, virialized filament at t = 73(kjc s ) . 



radial density p(r) and pressure P(r) profiles are given by 

LA. ( r l^L\ - _ 47rG 
r 2 dr\p dr ) " ' ' P cfl \p cfl 



\ d f r 2 dP\ A n P ( p \ 4 / 3 

4vrGp, — = [ — ) ■ (20) 
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Fig. 5. — Comparison of the density distribution for the conservative and non-conservative al- 
gorithm for simulations J and Jc in Table 1. The snapshots are taken at time t = 1.1, when 
the collapse stops and the energy error reaches its maximum value for the simulation using the 
non-conservative algorithm. The left panel is the density distribution from the non-conservative 
algorithm while the right panel is the result from the conservative algorithm. 

In our units, AttG = 1, the central density is p Cj o = 1, and the central pressure in hydrostatic 
equilibrium is P C; o = 8.09 x 10~ 5 . The adiabatic index 7 = 1.36 is chosen to be close to 4/3 in 
order to clearly demonstrate the difference between the two algorithms. The background density 
in the simulation is set to 0.005. The diameter of the sphere is 0.16, which is 20% the size of the 
simulation box, and the sphere is located at the center of the simulation box. The resolution we 
use is 256 3 . 

Initially, we decrease the central pressure to be a fraction of the pressure in hydrostatic equi- 
librium and let the sphere evolve under its own gravity. Let us estimate how much the size of 
the polytropic sphere should change due to a decrease in the central pressure. Starting from the 
virial theorem and using simple scaling relations, it is straightforward to derive that in hydrostatic 
equilibrium 

M 2_ 7jR 3 7 -4 = a ( 7 )P c0/9 -7 (21) 

Here a(j) is a constant that depends only on 7, R is the radius of the sphere, and M is the total 
mass. If we assume homologous adiabatic collapse, then the right hand side of equation (21) is a 
constant during the collapse. Suppose now that we reduce the initial pressure to be a fraction e of 
the pressure in hydrostatic equilibrium. We can then solve for the final radius, R, in terms of the 
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initial radius Rq as 



| = £ " < 22 > 



In our simulations, we use 7 = 1.36 and e = .93, which gives R/Rq ~ .4. Thus, we expect to get a 
large-amplitude, spherically-symmetric, periodic solution. Moreover, since the fractional radius of 
the sphere is small compared to the box size, we expect spherical symmetry to be maintained quite 
well, even though we use periodic boundary conditions. 

We start two simulations with exactly the same initial conditions but one uses the non- 
conservative algorithm whereas the other uses the conservative one. The time-evolution of the 
central density from the two simulations is shown in Figure 6. Both simulations give almost the 
same result for t < 15, but after the initial phase of the collapse, they show quite different behaviors 
with the non-conservative algorithm attaining a much larger peak central density. Snapshots of the 
density distribution from the two simulations through the plane y = at times t = 460, t = 680, and 
t = 900 are shown in Figure 7. Radial profiles of the spheres at time t = 460 and t = 900 through 
the line {x = 0, y = 0) are examined in Figure 8. At the early time, t = 460, both algorithms 
are able to keep the spherical symmetry of the object, although non-spherical structures have al- 
ready appeared with the non-conservative algorithm (top panels of Figure 7) . These non-spherical 
structures are amplified by time t = 680. Eventually, the initial density profile is destroyed for the 
non-conservative algorithm, and a low density hole is formed at the center of the sphere at time 
t = 900. On the other hand, the conservative algorithm (bottom panels of Figure 7) yields a peri- 
odic solution and keeps the spherical symmetry and initial profile of the sphere much better during 
the collapse than the nonconservative one. In conclusion, this test shows that for 7 ~ 4/3 in 3D, 
the dynamics are significantly affected by the choice of algorithm. The conservative algorithm is 
favored, since it better preserves spherical symmetry throughout the collapse and yields a periodic 
solution as expected. 



5. Performance 

In order to calculate the time derivative of the potential, we need to solve one extra Poisson 
equation compared to the old method. If the Poisson solver takes a significant fraction of the total 
simulation time, the code will be slowed down by almost a factor of two for this new feature. If 
the time spent in the Poisson solver is negligible, then the new code is almost as efficient as the old 
one. For all our simulations, we use periodic boundary conditions and solve the Poisson equation 
using FFTs. In this case, the new algorithm is almost as efficient as the original method. 
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Fig. 6. — Time evolution of the central density in the polytropic sphere test from two simulations 
that use the same initial conditions but different numerical algorithms. The solid line is the result 
using the non-conservative algorithm whereas the dashed line is the result using the conservative 
one. 

6. Summary and Discussion 

We have developed a new algorithm to conserve total energy to round-off error for Eulerian 
hydrodynamical simulations incorporating self-gravity. We have implemented this new algorithm 
in Athena and shown that it conserves total energy to round-off error as expected. By comparing 
a set of tests in ID, 2D and 3D with two different numerical algorithms for self-gravity (the new 
conservative algorithm and the traditional non-conservative algorithm), we have shown that the 
numerical error made in the traditional algorithm can change the dynamics significantly. From 
our numerical experiments, we conclude that the conservative algorithm will be important when a 
small amount of energy error can dramatically affect the dynamics, such as in radiation-dominated 
systems with an adiabatic index close to 4/3. 

We tested our new algorithm by implementing it in Athena, which uses an unsplit Godunov 
method to solve the equations of hydrodynamics. In principle, this algorithm can also be used in 
codes that implement an operator split scheme, such as ZEUS (Stone k, Norman 1992). In the 
application using ZEUS with self-gravity performed by Jiang & Goodman (2011) , the energy error 
using the traditional non-conservative algorithm was positive. If ZEUS is used to study problems 
of star formation with the non-conservative algorithm, gravitationally bound clumps can acquire 
a positive total energy due to numerical errors and be destroyed as a result of numerical heating. 
However, if our conservative algorithm is used, this will not happen, since the total energy is 
conserved. 

The new conservative algorithm is not necessary for all hydrodynamical simulations incorpo- 
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Fig. 7. — Snaphsots of the density distributions from the polytropic sphere test taken at different 
times through the plane y = from two simulations using different algorithms. From left to 
right, the times are t = 460, 680, 900, and in each case, only the central region is shown. The two 
simulations start from exactly the same initial conditions. The top panels correspond to the non- 
conservative algorithm whereas the bottom panels correspond to the conservative one. A central 
hole is created with the non-conservative algorithm and the simulation shows strong deviations 
from spherical symmetry. 

rating self-gravity. We find that for cases when the energy error is not a concern, the two algorithms 
give almost identical results. For example, for small oscillations of a polytropic sphere with an adi- 
abatic index of 7 = 5/3, the two algorithms yield very similar temporal and spatial behavior. In 
such cases, the old algorithm is just as good as the new one and is potentially more efficient. 
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Fig. 8. — Density profiles at different times on the line (x = 0, y = 0) from the 3D polytropic sphere 
test using two different numerical algorithms. The blue line is the initial density profile, which is the 
same for both simulations. The black lines are results for the simulation using the non-conservative 
algorithm whereas the red lines are results for the simulation using the conservative algorithm. The 
solid black and red lines are measured at time t = 460 while the dashed black and red lines are 
measured at time t = 900. 

A. Alternate forms 

If the only properties required of the gravitational energy flux are that it satisfy equation 
(11) and have the correct physical dimensions, then F g is not unique. Energy conservation is not 
affected by F g — > F g + a if V • a = 0. For example, 

a = — — Vx (6v x Vd>) . 
4ttG 

More generally, for arbitrary vector fields b, the replacements 

(F g ,E g )^(F g + d t b, E g -V-b) 

preserve dtE g + V • F g and therefore preserve the energy-conservation law. In particular, if one 
chooses b = — (167rG) _1 V(0 2 ), then starting from equation (13) for F g and \p§ for E g , one obtains 
the alternate forms ^ 

F 9 = p<Pv - , E g = pep + — \V4>\ 2 . (Al) 

Because they differ by a divergence, the two forms of E g yield the same total gravitational energy 
when integrated over all space. 

Additional criteria are clearly needed to select the "correct" forms of E g and F g . One such 
criterion is galilean covariance. The energy density, energy flux, momentum density, and stress 
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tensor of a nonselfgravitating fluid are 

E = ^pv 2 + u, F = (E+p)v, j = pv, T = pvv + pl, 

where u is internal energy per unit volume. Under an infinitesimal galilean transformation (v, p,p, u) — s> 
(v + SV, p,p, u), these densities and fluxes transform to first order in SV as 

E^E + j-SV, j^j + pSV, 

F^F + ESV + T-SV, T ->• T + 2j ■ SV . (A2) 

It seems reasonable to require that the gravitational contributions to E, F, and T should preserve 
the form of this transformation. Since the newtonian gravitational mass and momentum densities 
vanish under any sensible definition, E g and T g should be galilean invariants, while 

F g ^F g + E g SV + T g ■ SV. (A3) 

Because and <p are galilean-invariant, T g [equation (8)] and both forms of E g are invariant. 
However, <j) — >■ <f) — SV • V0, from which one can show that the forms (Al) of E g and F g are 
compatible with (A3), whereas the original forms E g = pcp/2 and F g as given by equation (13) are 
not. 

No matter what definitions are used for (E g , F g ,T g ), the same time evolution results if the 
conservative equations are solved exactly, provided that the definitions are mathematically equiva- 
lent to the original equations of motion (1) — (4). However, the spatial distribution of gravitational 
energy and stress will depend upon the inertial frame used, and it is possible that numerical trun- 
cation errors may be sensitive to this dependence. The whole point of the present exercise is recast 
the equations so as to improve energy conservation in finite-difference or finite- volume approxima- 
tions. Viewed in this light, numerical robustness becomes more important than formal elegance. 
We mistrust the form (Al) of the gravitational energy density because it is always non-negative 
in vacuum or near-vacuum regions: thus a small local error in the estimate of |V0| 2 /87rG could 
cause a large error in the update of the kinetic or internal energies per unit mass. The original 
form Eg = p4>/2 does not have this defect because it tends to zero smoothly with the mass density. 

If it is important both that the gravitational terms be locally galilean covariant and that 
E g = p4>/2, then we may keep equation (13) for F g but replace the gravitational stress tensor (8) 
with 

T a = ^G [( V< ^ V< ^ " ^ VV ^ + 2 P# • (A4) 
The difference between (8) and (A4) can be written as 

T a-^9 = t 1 V 2 " VV ] (A5) 

Since the divergence of (A5) vanishes, the gravitational momentum density remains zero if T g is 
replaced by T g . 
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